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ABSTRACT 

We present a new, simple, fast algorithm to numerically evolve disks of inelastically colliding particles 
surrounding a central star. Our algorithm adds negligible computational cost to the fastest existing 
collisionless N-body codes, and can be used to simulate, for the first time, the interaction of planets with 
disks over many viscous times. Though the algorithm is implemented in two dimensions — i.e., the motions 
of bodies need only be tracked in a plane it captures the behavior of fully three-dimensional disks in 
which collisions maintain inclinations that are comparable to random eccentricities. We subject the 
algorithm to a battery of tests for the case of an isolated, narrow, circular ring. Numerical simulations 
agree with analytic theory with regards to how particles' random velocities equilibrate; how the ring 
viscously spreads; and how energy dissipation, angular momentum transport, and material transport are 
connected. We derive and measure the critical value of the coefficient of restitution above which viscous 
stirring dominates inelastic damping and the particles' velocity dispersion runs away. 

Subject headings: accretion disks, planets:rings 
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1. INTRODUCTION 

How does a disk of collisional particles surrounding a 
star evolve in the presence of planets? The answer to this 
question has important implications. For example, after 
the planets of our Solar System accreted most of their 
mass, many small, rocky and icy bodies remained orbit- 
ing the Sun. Somehow, the planets eliminated most of 
these remnant planetesimals, while leaving some behind 
to form the asteroid belt, the Kuiper belt, and the Oort 
cloud. In the vicinity of Uranus and Neptune, the small 
bodies must have been highly collisional. Otherwise, these 
planets would have taken 10 12 yr to form in situ (Goldre- 
ich et al. 2004). 4 Yet virtually all simulations of the late 
stages of planet formation in the outer Solar System — such 
as those that model the migration of the ice giants, the 
resulting trapping of Kuiper belt objects into resonances, 
and the ejection of small bodies to the Oort cloud — neglect 
collisions. When the effects of collisions are accounted for, 
the current picture of the formation of planetary systems 
might change drastically. 

Planetary rings provide another setting in which inter- 
particle collisions play a crucial role. What are the origins 
of narrow rings shepherded by satellites? How do narrow 
rings settle into their special apse and node-aligned states 
(e.g., Chiang & Culter 2004)? And how do rings back- 
react upon and shape the orbits of shepherd satellites? 
Our understanding of satellite-ring interactions bears on 
mysteries such as the origin of eccentricities of extra-solar 
giant planets (e.g., Goldreich & Sari 2003). 

Despite its importance, the behavior of particle disks in 
the presence of perturbing bodies is poorly understood. 
Numerical simulations can help to further understanding. 
But until now, simulations of collisional disks have been 
too inefficient to follow, say, how disks viscously spread 
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in the long term. Collisions are traditionally simulated 
with a brute-force method (e.g., Brahic 1977; Wisdom & 
Tremaine 1988): at each time step of the integration of the 
gravitational equations of motion, it is determined which 
pairs of particles might collide before the next timestep. 
These potential collision pairs are then integrated forwards 
in time with a much smaller timestep, to see if they really 
do collide. But this method is inefficient: a brute-force 
search for collision partners requires around iV t 2 p opera- 
tions at each timestep, where N tp is the number of test par- 
ticles. In addition, most potentially colliding pairs do not 
collide, particularly in optically thin disks. Hence much 
computing time is wasted on missed collisions. More com- 
plex algorithms have been devised to reduce computing 
time (e.g., Lewis & Stewart 2000; Charnoz et al. 2001). 
But these are still not nearly as fast as the fastest colli- 
sionless N-body codes, such as SWIFT (Levison & Duncan 
1994). 

We sought a collision algorithm that (i) could be added 
to any N-body code, such as the freely-available SWIFT; 
(ii) contributes negligibly to the computational cost; (iii) 
is simple conceptually; (iv) is easy to code; and (v) follows 
correctly the long-term viscous evolution of disks in the 
presence of planets. We designed our algorithm to simu- 
late a vertically optically thin disk of identical, collisional, 
massless, inelastic but indestructible test particles that feel 
the gravity of the Sun and of multiple planets. Compli- 
cations that we do not include, such as the self-gravity 
of the particles, order-unity optical depths, and particles 
with differing sizes, spins, and cohesive strengths, could all 
affect the viscous evolution in ways that are not currently 
understood. But at this stage it seems wisest to ignore 
these complications, even though the algorithm could be 
modified to handle them. Viewed in its most basic terms, 
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inelastic collisions dampen random velocities and act as 
a source of friction between neighboring streamlines. As 
long as our algorithm preserves this behavior, while con- 
serving angular momentum and accounting for the loss 
of energy in inelastic collisions, it seems likely that it will 
properly model the long-term evolution of collisional disks. 
In the present paper, we test this assertion thoroughly 
when there are no planets, comparing in detail the results 
of our simulations with those of analytic theory. In a fu- 
ture paper, we shall test our algorithm in the presence of 
planets. 

2. THE COLLISION ALGORITHM 

The gravitational equations of motion of the Sun, plan- 
ets, and massless test particles are integrated with the 
Wisdom-Holman mapping method (Wisdom & Holman 
1991), using the SWIFT subroutine package (Levison & 
Duncan 1994). 

We supplement SWIFT with a subroutine that simulates 
collisions between test particles in a disk with vertical op- 
tical depth 
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that we actually use, since we apply this algorithm every 
time interval dt (and not i rb), we must correspondingly 
reduce the probability of a collision by di/t or b in order to 
maintain the collision time at the value given by Equation 
(2). 

Although carried out in only two dimensions, we empha- 
size that our algorithm models three-dimensional disks in 
which collisions maintain inclinations that are comparable 
to the random eccentricities. A truly two-dimensional disk 
is not realistic because collisions invariably generate out- 
of-plane velocities. But if one could somehow prevent the 
generation of out-of-plane velocities, the collision time in 
such a disk would be ~ s/(ut), which differs from Equa- 
tion (2) by the factor s/(u£ rb)- Since our algorithm satis- 
fies Equation (2), it does not model truly two-dimensional 
disks. 

As will be shown below, collisions drive the random 
speed of the particles toti> s gri d/iorb- Hence if two parti- 
cles fall in the same grid cell at one time, they will usually 
fall in separate grid cells after one orbital period. Since 
collisions can potentially occur every time step, it might 
be thought that the same two particles can collide many 
times in succession — a behavior that we consider undesir- 
able. But this behavior is avoided by the requirement that 
particles must be approaching each other for a collision to 
occur; immediately after they collide, their relative veloc- 
ity reverses sign, and they are no longer candidates for a 
collision pair. 

One of the main advantages of our algorithm is that the 
timestep is not restricted by the Courant condition. In a 
brute- force algorithm, one must restrict dt <C s/u in order 
to ensure that any two particles that fall within a dis- 
tance s of each other collide. This restriction on dt can be 
very cumbersome when u 3> s/i rb, as it will be whenever 
planets stir up the eccentricities. We avoid the Courant 
condition by treating the vertical dimension statistically: 
when two particles fall within the same two-dimensional 
grid cell, they need only collide a small fraction of the 
time because their vertical positions will, in general, dif- 
fer. With our algorithm, wc may choose dt to be as large 
as is allowed by SWIFT, which is typically a significant 
fraction of the orbital time. 

If two particles have been selected for a collision, i.e., if 
they lie in the same grid cell, are approaching each other, 
and are selected by the random number generator, then 
their velocities are updated as though the bodies were fric- 
tionless spheres whose surfaces touch (e.g., Trulsen 1971): 
the component of the relative velocity vector that lies par- 
allel to the axis connecting the two particles is reversed in 
sign (from a converging velocity to a diverging one), and 
multiplied by the coefficient of restitution e, i.e., in obvious 
notation, 

<4el,|| = - eU rcl,|| • (3) 

Neither the perpendicular component of the relative ve- 
locity vector, nor the velocity of the center of mass of the 
two colliders, nor the positions of the colliders are changed 
by the collision. A collision does not alter the sum of the 
angular momenta of both colliding bodies; hence the colli- 

5 More precisely, in optically thin disks the r.m.s. azimuthal speed is twice the r.m.s. radial speed; the r.m.s. vertical speed is comparable. 

6 Although this simpler algorithm yields the correct collision time, we did not use it because it introduces an artificial frequency into the 
problem, set by the time interval at which the algorithm is applied (~ t or b)- When wc attempted this algorithm, wc found that a gap was 
cleared in the disk of test particles where the orbital period was exactly equal to this interval. 



where N tp is the number of test particles, s is their size, 
and f and A are, respectively, the mean orbital radius 
and the radial width of the annulus that the particles 
occupy. In collisional particle disks, collisions tend to 
isotropize the velocity distribution. 5 The collision time 
is t co \ ~ l/(n v s 2 u), where u is the 1-D random speed and 
n v is the volumetric number density, which is related to r 
via r <~ n v s 2 ut OI b. Therefore 
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The collision time is longer than the orbital time by the 
u- independent factor 1/r. 

We capture this behavior with two-dimensional simula- 
tions in which all bodies have zero inclination. Every time 
step dt, a two-dimensional square grid is built, with each 
grid element having dimensions s gr id x s gr ;d; s gr id can be 
thought of as the size of a particle. If two test particles fall 
in the same grid cell, and if their relative speed is negative 
(i.e., if they are approaching each other), then they col- 
lide with each other with probability P co i = dt/t OT b <C 1, 
where i or b is the orbital time at the collision point. A ran- 
dom number generator is used to determine whether or 
not they actually collide. 

To sec that this algorithm gives the same collision 
time as Equation (2) (where t is given by Eq. [1] with 
s — > s gr id), it is instructive to consider first a simpler al- 
gorithm that also yields the correct collision time. In this 
simpler algorithm, one waits for a time interval of t or b (in- 
stead of dt) before finding which particles fall in the same 
grid cell. Then two particles which do fall in the same grid 
cell, and have converging velocities, collide with probabil- 
ity -Pcoi = 1- Since the probability that a given particle lies 
in a cell occupied by a second particle is r, the collision 
time is t OT b/r, as required. 6 Turning now to the algorithm 
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sion algorithm exactly conserves total angular momentum. 
Note that a collision between two particles separated by 
distance d changes the velocities of the particles as though 
each was a smooth sphere with radius d/2. Since d changes 
from collision to collision, the particles' sizes arc effectively 
changing; they are only approximately s gr ;d- 

The algorithm has now been completely described, aside 
from how the code finds which pairs of particles lie in the 
same grid cell. To find colliding pairs, the code first de- 
termines in which grid cell each particle lies. A grid cell 
is labelled by two integers, representing its location along 
the x- and y-axes. Second, the code sorts the grid cells 
that contain test particles with the heapsort algorithm 
(Press ct al. 1992). The sorted occupied grid cells are then 
checked to see if the same grid cell is repeated for two dif- 
ferent particles. The step that takes the most time in the 
entire collision algorithm is the heapsort, which requires 
~ iVtplnA^tp operations. However, in the runs presented 
in this paper, with N tp — 10 4 particles, it was found that 
the collision algorithm contributed negligibly to the run- 
ning time of SWIFT. 

3. SIMULATIONS OF NARROW CIRCULAR RINGS 

In the remainder of this paper, we investigate circular 
rings of particles, without any planets. Circular rings are 
understood quite well theoretically (e.g., Lynden-Bell & 
Pringle 1974; Brahic 1977; Goldreich & Tremaine 1978; 
Shukhman 1984; Petit & Henon 1987). Our goal is not 
only to test the collision algorithm, but also to develop di- 
agnostics that can be used for the much more complicated 
case when planets are present. 

The parameters of the simulation include the coefficient 
of restitution (e), the size of a grid element (s gr ;d), the 
number of test particles (N tp ), and the initial orbital el- 
ements of the test particles. The central body's mass is 
IMq. We simulate narrow rings with mean radius f = 1 
AU and radial width A <C f, and choose dt = 0.18 years 
for the integration timestep. 

A particle ring has two characteristic timescales (Brahic 
1977). The shorter one is t co \, the time for a particle to 
collide. The longer one is t^is, the time for particles to 
diffuse across the ring's width. Random velocities relax to 
their equilibrium distribution on timescale t co \, while the 
ring density evolves on timescale idiff- We investigate in 
turn the evolution of random velocities and of density. 

4. RANDOM VELOCITY EVOLUTION 

4.1. Theory 

Collisions can both excite random velocities by draw- 
ing energy from the background Keplerian shear ( "viscous 
stirring" ) , and damp random velocities because of finite in- 
elasticity. When the coefficient of restitution e < e*, where 
e* is a critical value, a typical collision damps random ve- 
locities. But there is a limit to how cold the particles get. 
Particles on circular orbits collide with one another at the 
Keplerian shearing speed ~ sfl, where s is the size of a 
particle and O is the orbital angular frequency. In our col- 
lision routine, this relative speed is ~ s gr jd^- Since a single 
such collision redirects the particles onto non-circular or- 
bits, the random velocity cannot fall below ~ s gr idf2, and 
the rms (root-mean-squared) eccentricity always relaxes to 

e rms ~ ^ (4) 



in the e -C e* limit, where r is the local disk radius. 

By contrast, when e > e_, a typical collision excites ran- 
dom velocities: viscous stirring dominates inelastic damp- 
ing, and the rms eccentricity runs away. 

In optically thin disks composed of equal-size parti- 
cles, e* is determined solely by the angular dependence 
of the differential collisional cross-section. For frictionless 
spheres, e* = 0.63 (Goldreich & Tremaine 1978). For fric- 
tional spinning spheres, e* = 0.92 (Shukhman 1984). For 
our collisional cross-section, we derive e* = \/7/5 = 0.529 
(Eq. [69]). 

4.2. Simulations: Approach to Velocity Equilibrium 

Figure 1 shows the time evolution of the rms eccentric- 
ity 

e rms = (e 2 ) 1 / 2 , (5) 

in three simulations, where () averages over all test parti- 
cles. In each simulation, s gr id = 10 -3 AU and N tp — 10 4 . 
Particles were initially on orbits evenly spaced in semi- 
major axis between f — Ao/2 and f + Ao/2, where f = 1 
AU and A = 0.08 AU (top hat profile). Initial eccentric- 
ities were identical, and initial longitudes and pericentre 
longitudes were random. 
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Random Velocity Evolution. The two simulations with 
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e = 0.3 both relax to the Same 6rms ^ 
simulation, heating proceeds indefinitely. 

Two simulations had e = 0.3. One of these was initially 
cold, with initial eccentricities = 0; the other was initially 
hot, with eccentricities = 0.01. In both simulations, e rms 
approached <~ s^aff = 10~ 3 (Eq. [4]). Clearly, inelastic 
damping dominates viscous stirring when e = 0.3. The 
third simulation had e = 0.6, and was initially cold. Its 
e rms grew indefinitely. Hence viscous stirring dominates 
when e = 0.6. 

We define the collision time as 
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In the three simulations, we measured t co i = 106 — 108 
yr. An estimator of t co \ from the input parameters of a 
simulation follows from a more precise form of Equations 
(1)"(2): 



where 
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is the number of particles per radial distance, and 4w is 
the product of two factors: 2n for the area of a ringlet 
(2irr x dr), and 2 because only half of the time, when the 
relative velocity is negative, do doubly-occupied grid cells 
lead to collisions. For the parameters of the present simu- 
lations, with r = f and dN tp /dr — N tp /A , = 101 yr. 

Because t^oi*' ^ oes n0 ^ accoun t f° r the decreased collision 
frequency of particles at the edge of the ring, it underesti- 
mates t co i by a small amount. That t co \ remained nearly 
constant throughout the simulation reflects the near con- 
stancy of the ring's width, since parameters were deliber- 
ately chosen to freeze out diffusion (from Eq. [9] below, 
the diffusion time in these simulations is ~ 10 5 yr). 

In the cold e = 0.3 simulation, e rms reached 80 percent 
of its final value by time t — 1000 yr — around 10 colli- 
sion times. The hot e = 0.3 simulation took much longer 
to reach velocity equilibrium: initially, e rms decayed ap- 
proximately exponentially with a time constant of 2500 
yr — around 25 collision times. Even after t = 2500 yr, the 
hot simulation took hundreds of collision times to reach 
its final e rms . Velocity equilibration in the hot simula- 
tion was so long because of particles at ring edges. Edge 
particles tend to retain their initial eccentricities because 
their epicyclic excursions carry them away from the ma- 
jority of particles; consequently, edge particles collide less 
frequently than do particles in the ring proper. 



Fig. 2. — Equilibrium e rms . Seven points denote the steady- 
state e rms in simulations with differing e. Vertical line shows crit- 
ical value e* = V7/5 = 0.529. In initially cold simulations with 
e = {0,0.1,0.2,0.3}, e rms was evaluated at 10 4 yr. In initially cold 
simulations with e = {0.4, 0.5, 0.525}, e rms was evaluated at 2 X 10 4 
yr, since these simulations took longer to reach steady state. 

4.3. Simulations: Equilibrium Eccentricities 

Figure 2 shows results from seven simulations, all with 
the same initial conditions as the cold simulations de- 
scribed in the previous subsection, except with differ- 
ing values of e. The seven plotted simulations, with 
e < 0.525, all reached velocity equilibrium, with 
(order unity constant) x s gr id/^- Simulations with e > 0.6 
never reached velocity equilibrium. We conclude that 
0.525 < e* < 0.6 for our collisional cross-section. In Equa- 
tion (69) below, we derive e» = V7/5 = 0.529. 

5. DENSITY EVOLUTION 

5.1. Theory 

A ring diffuses in the time that it takes a particle to 
random walk across its width. This random walk has 
a step-size equal to the epicyclic excursion of a particle 
(~ re rms ~ Sg r id) arid a time per random step of t co \. So 
to diffuse the width of the ring A takes a time 



tdiS ~ tcoi 



A 



(9) 



where the inequality holds when A » s gr id; otherwise, 
t<us ~ tco\ until A <~ Sg r id- Since t col cx 1/n oc A (Eq. [7]), 
a ring expands as 

Aoct 1 / 3 . (10) 



More precisely, n satisfies the diffusion equation 
dn d ( dn 



dt dr \ d 



(11) 



where the viscosity v — constant x Sg rid /t co i (Eq. [9]). 
Inserting Equation (7) into our expression for v, we see 
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that 



v = k v 



grid 
^orb 



n 



(12) 



which defines the dimcnsionlcss constant fc„; k v is a func- 
tion of e, but is independent of s gr id, r, and n. Pe- 
tit & Henon (1987) derived the above diffusion equation 
and gave its self-similar solution, an expanding inverted 
parabola: 



n = 



where 



2 A 




\r-f\< A/2 (13) 
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A = 36k„ 



grid 
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Since the viscosity v decreases with decreasing n, the diffu- 
sion is non- linear, and the edges of the ring at r = f ± A/2 
are sharp. 

5.2. Simulations 

Because of the steep dependence of the diffusion 
timcscalc on the width of the ring, t oc A 3 , it takes a 
long time to simulate even a modest increase in A. Sim- 
ulation parameters must be chosen judiciously. We fix 
N, 



Up 

for s 



10 4 and r = 1 AU, and seek the optimal values 
Kri d and Aq = A|t=o. To simulate as large an in- 



crease in A as possible, the simulation should begin with 
as narrow a ring as possible. For fixed s gr id> the narrow- 
est ring that is not optically thick has unity optical depth: 
A <~ A r tp(s gr id) 2 /27rf. The evolution timescale at the 
start of the simulation is to = constant x (A ) 3 /(s g rid) 4 
(Eq. [14]), so with unity optical depth, t = constant 
x Ao. Hence the fastest timescale is obtained with the 
smallest Ao. But we must have Ao > s g rid 7 so the optimal 
values are Ao = s gr id = 2irf/N tp . Rounding up, we set 

A = s gr id = 10_3 AU - 

Figures 3 and 4 portray results from a simulation with 
the parameters listed above, and e = 0.3. Initially, the ring 
particles were uniformly distributed in a ring with edges at 
1 ± (5 x 10~ 4 ) AU. In Figure 3, we show how the particles' 
dispersion in r, 



((r-m 



2\l/2 



(15) 



and how the collision time (t co \) vary with time, according 
to both theory and simulation. Figure 4 displays compar- 
isons between theory and simulation for n. According to 
the theory described in §5.1, n(r) is given by Equation 
(13), the dispersion in r is 

o-l = j(n/N tp )(r-f) 2 dr = A 2 /20 (16) 

with A given by Equation (14), and the collision time is 



1207rife„^gf£i . (17) 



The theory fits the numerical results well when we set 

k v = 0.016. (18) 
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Fig. 3. — Dispersion and Collision Rate Evolution in a Diffusing 
Ring. Top panel: points show oy = ((r — f) 2 } 1 / 2 from the numerical 
simulation. The theory line through the points is A/\/20 (Eq. [16]), 
where A is given by Equation (14). The normalization of the line 
was adjusted by choosing k v = 0.016. Bottom panel: points show 
the collision time t co \ from the simulation (Eq. [6]). The theory line 
is given by Equation (17), with fc„ = 0.016. Theory underestimates 
*col by a small amount because of particles at ring edges (§4.2). 

100 



80 



h 60 



s 



40 



20 



1 1 1 1 1 1 


■ i 1 


I I | I I I | I 

e = 0.3 


t = 10 3 yr 










[ t = 10 s yr '_ 


" , A, , , , 1 


, 1 , 





0.96 



0.98 



1.00 

r(AU) 



1.02 



1.04 



Fig. 4. — Density Evolution in a Diffusing Ring. Histograms of 
the number density are shown at two times from the same simula- 
tion as in Figure 3. The theory lines through the histograms are 
given by Equation (13), with fc„ the same as in Figure 3. 

6. ANGULAR MOMENTUM AND ENERGY TRANSPORT 

Understanding how planets interact with disks requires 
understanding how angular momentum and energy are 
transported — both within disks and between planets and 
disks. Below we study transport in isolated, circular, nar- 
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row rings, developing diagnostics that will prove useful in 
future simulations of disks with planets. 

6.1. Theory 

The azimuthally-averaged equations describing the con- 
servation of particle number, angular momentum, and en- 
ergy are 

8 t n + r F n = (19) 

d t (Hn) + d r (HF n ) = -d r Fji sc (20) 

d t (En) + d r (EF n ) = -8 r F E isc - nt , (21) 

where F n is the (net) number flux across a circle of radius 
r, i.e., the number of particles per unit time that exit this 
circle minus the number that enter it; H(r) = (GMqt) 1 / 2 
is the specific angular momentum of a particle on a circu- 
lar orbit; E(r) = -GM Q /2r is the specific energy; FJj sc 
(the "viscous flux of angular momentum" ) is defined as the 
difference between the total angular momentum flux and 
HF n , i.e., Fp sc 

difference in energy fluxes, i.e., F E 1SC + EF n 
is the rate at which specific energy is lost, per particle, in 
inelastic collisions. 

Our decomposition of F^ into two components has the 
following interpretation: HF n is the angular momentum 
flux that would be carried by particles on circular orbits 
whose radii change on timescales long compared to the or- 
bit time, while FJj sc is the part of the angular momentum 
flux not associated with the direct advection of circular 
orbits. The viscous flux of angular momentum and the 
corresponding viscous flux of energy are transferred in the 
ratio appropriate for circular orbits, i.e., since circular or- 
bits have E = -( 1/2) (GM Q /iJ) 2 , a transfer in angular 
momentum of 5H must be accompanied by a transfer in 
energy of SE = {GM Q ) 2 SH/H 3 ^il-SH, so 



„ + HF n = F^ ot ; F| isc is the corresponding 

F| ot ;and£ 



jpvisc 
E 

TTivisC 

H 



= Q(r) . 



(22) 



The above relation applies only to the viscous fluxes, not 
to the total fluxes (i.e., EF n /HF n = — £1/2). Equations 
(20) and (21) simplify with the aid of Equations (19) and 
(22) to 



(23) 

F£ sc = -n£/(dQ/dr) . (24) 

The latter is the well-known relation between energy dis- 
sipation and viscous angular momentum flux for accretion 
disks (e.g., Lyndcn-Bcll & Pringle 1974). 

Since £ ~ s 2 rid Jl 2 /i co i ~ ns* lid fl 2 /4nrt olh , we may re- 
write (19) using (23) and (24) as 



where 



d t n 



2rt n 



,2 a 4 



S «rid( 3fi2 /8rtorb) 



(25) 



(26) 



is an e-dependent dimcnsionless constant. In deriving (25), 
we have dropped terms that are small for narrow rings 
(e.g., \d(lnfl)/dr\ < \d(lnn 2 )/dr)\). Equation (25) is iden- 
tical to Equations (11)— (12) provided k E = k v . 



6.2. Simulations 

In this subsection, we diagnose angular momentum and 
energy transport in numerical simulations and compare to 
theory. To measure the number flux F n across a circle of 
radius r, we evaluate N >r , the number of particles whose 
radial distances from the Sun exceed r, at two times, t m 
and t m + dt m . Then 



F n {r : t m ) 



N >r (t m + dtm) - N >r (t m ) 



(27) 



where the symbol = means that this is how F n is mea- 
sured. Similarly, we measure the total angular momentum 
flux as 

H>r(tm + dtm) ~ H >r (t m ) 



F%\r,t m ) = 



dt r 



(28) 



where H >r is the sum of the specific angular momenta 
of all particles whose radial distances exceed r. For the 
energy flux, we must account for the energy lost in inelas- 
tic collisions (where the specific energy lost per collision 
= m 2 c1 || (1 — e 2 )/4; sec Eq. [3]). To this end, when calcu- 
lating the energy flux across a circle with radius r, we first 
calculate 5£ >ri the total specific energy lost between times 
t m and t m + dt m in all inelastic collisions that occurred at 
radii > r. Then the total energy flux is 

E> r (t m + dt m ) — E >r (t m ) + S£>r 



Hr t 



dtr, 



(29) 



where E >r is the sum of the specific energies of particles 
with radii greater than r. The viscous fluxes are deter- 
mined by the three fluxes above: 

HF n (30) 

EF n . (31) 

Energy dissipation in an annulus between radius r\ and ri 
is measured via 



FiVISC 
H 
rpvisc 
E 



rptot 
E 



n£ = 



(32) 



(r-2 - n)dU 

Figure 5 shows a number of measurements of the angu- 
lar momentum flux for the simulation whose parameters 
are given in §5.2. The theory curve is from Equations (24) 
and (26): 



FiVISC 
H 



k E 



,2 4 , 
grid 



(33) 



with n given by Equation (13), and k E — k v = 0.016 (we 
verify the equality of k E and k v in Figure 7 below). Over- 
laid on this theory curve are three sets of datapoints, mea- 
sured in three independent ways. The agreement between 
theory and simulations is excellent. 

Figure 6 shows the number fluxes for the same simula- 
tion. The theory line shows (Eqs. [23], [33]) 



F n = -k t 



b grid 



On 2 



2rt orb dr 



(34) 



with n from Equation (13). The datapoints agree well with 
the theory, although there is some scatter. 

The circles in Figure 7 show the energy dissipation con- 
stant k E (Eq. [26]) for the simulations described in §§4.2- 
4.3, at the same output times (Fig. 2). Recall that each of 
these simulations has a top hat density profile and is run 
for much less than a viscous time, so the density hardly 
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evolves. The simulation with e = 0.3 has &e = 0.016. 
Therefore ks = k v (Eq. [18]), as suggested below Equa- 
tion (26). Note that it is much more efficient to measure 
kE than k v , since ks can be measured in only a few col- 
lision times, whereas k v must be measured on the viscous 
timescale. We defer to §6.4 a discussion of the diamonds 
and stars in Figure 7. 



< 



£ 1 



1 1 1 1 , 1 


i i i , i i i 1 


■ • F™ c 


□ -n<S/(d£2/dr) - 




— F™ c theory 




t m =l.lxl0 5 yr : 


, (?) , , 1 , 


■ ■■■■■■ .Hvi , 



0.96 0.98 



1.00 

r(AU) 



1.02 



1.04 



Fig. 5. — Angular Momentum Flux Measured With Various 
Methods. Data are taken from the simulation described in §5.2 
(Figs. 3-4) at time t m = 1.1 X 10 5 yr, with measurement interval 
dt m = 10 4 yr. The theory line is Equation (33). FJj BC is measured 
with Equation (30). Fg 1BC is measured with Equation (31); clearly 
F E ' BC /Q = FJj BC , confirming Equation (22). nS is measured with 
Equation (32); the squares confirm Equation (24). 



U 



fa 




1.00 

r(AU) 

Fig. 6. — Number Flux. Data are taken from the same simulation 
as in Figure 5, and at same times. The theory line is Equation (34). 
Datapoints are measured with Equation (27). 
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Fig. 7. — Energy Dissipation. The simulations are the same as 
those described in Figure 2. Circles show kE, measured with Equa- 
tions (26) and (32). Measurement parameters include n = 0.99 
AU, r 2 = 1.01 AU, dt m = 5000 yr, and t m = 5000 yr for 
e = {0, 0.1, 0.2, 0.3} and t m = 15000 yr for e = {0.4, 0.5, 0.525}. For 
e = 0.3, fc_E = 0.016, a value that matches k v as given by Equation 
(18). Diamonds show k r ^ <x F r( p, measured as described in §6.4. 
Stars give k 



NL 

E 



kE — k T Jt oc Fnl; this quantity is practically con- 



stant with e, as expected from Equation (52). Most of the energy 
dissipation (viscous transport of angular momentum) arises from 
F r and not from Fnl as e approaches €*; compare with Equations 
(51)-(52). 

6.3. Dynamics from a Microscopic Perspective: Theory 

The angular momentum flux advected by particles is 

= n(r(Clr + v^Vr) (35) 

= HF n + nr(v r v,p) , (36) 

where v — (v r ,v^,) is the difference between a particle's 
total velocity (in the radial r and azimuthal <j) directions) 
and the Keplerian circular velocity at its position, and () 
denotes an average over particles in a narrow ring. Defin- 
ing 

F r(l , = nr(v r V4,} , (37) 

we have 

Fjj 1 = HF n + F r< j, + -Fnl , (38) 

where _Fnl is the "non-local" flux, i.e, the angular momen- 
tum flux not advected by particles (Wisdom & Trcmaine 
1988). In §6.1, we considered only the combination 

^ isc = F r<p + F NL . (39) 

In this subsection, we wish to calculate F™ c in terms of 
microscopic quantities. Hence we must consider the two 
components of F™ c separately. The two components have 
the following interpretation: (i) A particle that crosses a 
circle of radius r has an angular momentum H' that is 
not, in general, equal to H(r); the difference H 1 — H(r) 
contributes to F r< p. (ii) When a particle that is inside of 
the circle collides with one that is outside, the angular mo- 
mentum transferred across the circle contributes to Fnl- 
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For a collisionless Kcplerian particle, to lowest order in 



e, 



v r = erfl sin(ftt) 
«0 = ^erft cos(flt) 



(40) 
(41) 



with r, il independent of time. Therefore for collisionless 
particles, F r( f, cx (w r «0) cx (sin(2f2t)) = 0. But collisions 
give a definite contribution to (vrV^) (Greenberg 1988), as 
we presently show. 

Instead of averaging over space only (as denoted by 
(VrVt/,)), it will prove convenient to average over both space 
and time. We average over a radial width Ar that is larger 
than the particle size but smaller than the lengthscale over 
which the density varies. We also average over a time At 
that is longer than a few collision times but shorter than 
the timescale over which the density changes. The average 
F r< j, is defined via 

/t+At ,-r+Ar 
J dt'drinr'vrv^) , (42) 

where v is a Lagrangian quantity that is tied to particles. 
Extracting F rt j, from a numerical simulation with only a 
spatial average ((iyt^)) can lead to large errors. In par- 
ticular, in an optically thin disk, only a small fraction of 
the particles have collided within the last orbital period. 
Hence most particles contribute negligibly to {v r v^,}, and a 
calculation of {v r v^} can become plagued by small-number 
statistics. 

For a single particle that does not collide in the interval 
At, 



t+At 



v r v<pdt' 



-ft 



J dt 



But if it collided once, 

r t+At 



-n 1 (vlu+At -v%\t) 

(43) 



dt' 



-O- 1 (v 2 



I t+At 



vl\t 



S(vD) , (44) 



where <5 represents the change due to the collision. There- 
fore each collision contributes to Equation (42) in the 
amount of 

(AiAr) 5F^ = rQ" 1 ^,!) + 5« 2 )) , (45) 

where 5(%n -J, 5{tn 2 ) are the contributions from the two 
collision partners. The contribution from the endpoints, 
v| I t+At — v^\t, can be neglected as long as At is much 
longer than the collision time. 
We define 



(AtAr)F 



NL 



dt' dr' _Fnl 



(46) 



If two particles collide when their positions are at radii 
r\ , r 2 (where r < r 2 < r± < r + Ar) , and if the particle 
at n has its </>- velocity changed by 5v$ t i in the collision, 
then the collision contributes to the integral in Equation 
(46) in the amount of 

(AtAr) <£F NL = (n - r 2 )r 1 6v <t , tl (47) 
wr(n -r2)(^,i-H,2)/2. (48) 
To obtain the latter symmetric form, we approximated 
ri ~ r and 5v<j, t i « —Sv^ >2 (when in fact r-\_5v^, t \ = 
—r 2 Sv < p. 2 )- The error accrued is of order (n — r 2 )/r ~ 
s/r -C 1, where s is the particle size. 



Since the number of collisions per unit time per unit 
radius is n/2t co \, 



2r 



r<fc 



-Fnl — 



*-col 

n r 



(49) 

, . . , ...„. (50) 

^col ^ 

where (} c (not to be confused with ()) is an average over 



((n - r 2 )5(v<f, t i - v^ 2 )) c 



collisions, and (5(v?)} c 



2 )) c /2. With our 



collision algorithm, t co \ may be pulled out of_the averages. 
We use Equations (49) and (50) to extract F rt p and Fnl 
from the simulations (see Figure 8 below). 

We can estimate the magnitudes of the two fluxes 
as follows. Since the peculiar velocity distribution is 
anisotropic, with = (v 2 )/4 < (v 2 ) (Eqs. [40]- 

[41]), and since collisions tend to isotropize the distri- 
bution, therefore collisions systematically increase v 2 , by 



(S(vl)) c 



+e 2 ms (rf2) 2 , transporting of F r< p out- 



wards. The contribution to Fnl is ((ri — r 2 )d(v<f,^i — 
v 4>,2))c ~ +fis 2 1 where s is the particle size. This is true 
even when e rms 3> s/r (i.e., when e — » e»), because in 
that case S(v<f, t i — v^ t2 ) is nearly random and hence nearly 
uncorrelated with n — r 2 . But because of the mean Ke- 
plerian shear, a small correlation ~ Vis 2 persists. This 
contribution also transports angular momentum outward. 
In sum, 



■ rfl^—(re r 

tco\ 

rfl—s 2 . 



)2 



(51) 
(52) 



6.3.1. Relating n£ to F% sc 

We now calculate the energy lost in an inelastic colli- 
sion, and thereby re-derive Equation (24) from a micro- 
scopic perspective. Consider the collision of two parti- 
cles having total velocities V \. 2 = V\, 2 + V C " 2 C , where 

= ^(ri, 2 )ri i2 4> is the circular Keplerian speed at the 
positions of the particles. Then the specific energy lost per 
collision is, in previous notation, 

(AtAr) 5 = -S(V? + V 2 2 )/2 (53) 

= -S(v 2 + v 2 )/2 

-S Vl -Vf c - Sv 2 -Vf c (54) 
= -S(v 2 1 +v 2 2 )/2 

-{r\ — r 2 )ri8v ( j >t idQ./dr (55) 
«-<*(«? + «|)/2 

-(n - r 2 )r(dtt/dr) x 
{5vt,i - Sv^/2 , (56) 
where to derive (55) we used conservation of orbital angu- 
lar momentum in a collision {r\5v$ t \ = — r 2 <5u^ 2 ). There- 
fore 



nE 



,, 2 s >' dfl . 



> /: "M"*) + 2 ^r( r i - r 2)^Ki - «0, 2 )>c • (57) 

We re- write the first term using v 2 = — 3w^ + e 2 (rJ7) 2 (Eqs. 
[40]-[41]). Since the particles are in collisional equilibrium, 
(% 2 )>c=0, (58) 
which implies that (S(v 2 )) c = —3(S(v^)) c . Inserting this 
into Equation (57), and comparing the result with Equa- 
tions (49)-(50), completes our proof of Equation (24). 
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6.3.2. Calculating F r< j, and e* 

We calculate the numerical constant that we dropped 
in Equation (51), and then use that result to calculate 
e*. Goldreich & Trcmaine (1978) perform a similar cal- 
culation, with a different (though still idealized) colli- 
sional cross-section. Shukhman (1984) accounts for Fnl as 
well. The treatments of Goldreich & Trcmaine (1978) and 
Shukhman (1984) arc more rigorous than ours, as they in- 
tegrate over the velocity distribution function. They also 
consider the more general case of disks with order-unity 
optical depth. But their final expressions are "extremely 
cumbersome" (Shukhman 1984). Although our treatment 
is not rigorous, it is considerably simpler, and we justify 
it by comparing with simulations. We neglect Fnl in the 
present subsection, taking e rms 3> s/r. 

For F r< p (Eq. [49]), we seek 

(S(vl)) c = (vl) ac (vl) bc , (59) 

where ()f, c is an average over collisions of the state imme- 
diately before the collision, and () ac is of the state imme- 
diately after. In an optically thin disk, the averages over 
particles arc equal to time-averages for a single particle: 



= -^i(rft) 2 , (vl) = 



(60) 



using Equations (40)-(41). With our collision algorithm, 
the probability that a particle collides is uniform in time. 
Therefore 

( V {r,<t>})bc = ("{r,*}) • ( 61 ) 

In collisional equilibrium (Eq. [58]), 

(S(v 2 r +4v 2 )) c = . (62) 

We now evaluate (v 2 ) ac and {v^) ac . We make the plausible 
assumption that the relative velocity of collision partners, 7 
u = vi v 2 , is isotropically distributed after the collision: 

{u 2 r ) ac = (u 2 ) ac . (63) 

This assumption is verified by numerical simulation in 
§6.4. To relate u to v, 

1 \c = 2(v\ rM ) bc (64) 

(65) 



W«{r,*})>c 



For the first relation, we neglected the correlation between 
Vi and V2 before a collision, and for the second rela- 
tion, we used S(vi + v 2 ) = 0. Equations (63)-(65) yield 
(S(vl - v 2 )) c = (l/2)(v 2 - vfybc, which, with Equations 
(60)-(62) becomes 

<<H^)>c = J^ 2 ms (rf]) 2 , (66) 



giving the numerical constant for F rt j, (Eq. [49]). 

To calculate , we evaluate the specific energy lost per 
collision in the ccntcr-of-mass frame: 



\{5{u 2 ) 



= A e 2 ms(r fi)2 

80 rmsV ; 



(67) 



(Eqs [62], [65], [66]). We can also evaluate the energy loss 
as follows: from Equation (3), it is 



e2 )< u rel,||/6c 



) bc = (l-e 2 ) 1 -(u 2 ) 



be 



(68) 



where we have made use of the isotropy of the pre-collision 
u with respect to the axis connecting the centers of the two 
particles. Since (u 2 )bc = (5/4)e 2 ms (rf2) 2 , Equations (67) 
and (68) are equal only if e is equal to 



^ = 0.529 

5 



(69) 



This is interpreted as the one value for e that enables 
the velocity distribution to equilibrate in the limit that 
> s/r, i.e., F r< p > _F NL . 



6.4. Dynamics from a Microscopic Perspective: 
Simulations 

In Figure 8 we test the theoretical results derived in 
§6.3. In the top panel, for the simulation with the same 
initial conditions as the one with e = 0.525 in Figures 
2 and 7, we compute (S(v,p) 2 ) c — the key factor that en- 
ters into F r( j, (Equation [49]) — in two ways: first as an 
average over all particles that collide within successive 
time intervals of 2000 yr, and second using Equation 
(66), where e 2 ms = (e 2 ) is a spatial, not temporal, aver- 
age over all particles (regardless of whether they collide). 
The agreement verifies Equation (66). Also shown in the 
top panel is the corresponding contribution to Fnl, i.e., 
(fi/4)((n -r 2 )S(v 4>>1 -tty, 2 ))c (Equation [50]). That this 
quantity is constant with time and sits far below (c>(i\/,) 2 ) c 
is expected from Equations (51)— (52), given this simula- 
tion in which e rms > s glid /f (i.e., e* — e <C e*). 

In the bottom panel of Figure 8, we test our assumption, 
made in Equation (63), that post-collision relative veloc- 
ities are isotropic. The measured near-equality between 
(u 2 )ac an d {vfyac ^ s satisfactory. 

Finally, returning to Figure 7, we plot separately the 
contributions to the total energy dissipation from F r ^ and 
i*NL- The former contribution is described by 



VT<P — 
K E — 



-F r4> (dn/dr) 



n 2 4 id (30 2 /8rt 

orb ) 



(see Equations [24], [26], and [39]). We measure F r , 
cording to 



>ri 



dtmin - r 2 ) 



(70) 



ac- 



(71) 



where 5(v^) >r is the total change in v 2 ^ summed over par- 
ticles that collide between times t m and t m + dt m , at radii 
> r. For simplicity, we evaluate fc^ L = ki 



L.r<p 



7 Since we take e > s/r, we neglect the difference in the Keplerian circular velocities at the positions of the two particles. 
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Fig. 8. — Averages over Collisions. The two panels show the time 
evolution of various collisional averages in a simulation with the 
same initial conditions as the one with e = 0.525 in Figures 2 and 7. 
Top panel: open diamonds show the factor that enters into F ri p (Eq. 
[49]), i.e., {S(v^)) c = (l/2)(<5(ir| j + v 2 ^ 2 )) c , normalized as shown. 
The averaging over collisions is done by recording the peculiar veloc- 
ities of colliding particles immediately before and after each collision. 
The averaging time is 2000 yr. Solid diamonds give e 2 ms = (e 2 ) (a 
spatial, not temporal, average over all particles), multiplied by the 
appropriate prc-factor as given by Equation (66). The agreement be- 
tween open and closed diamonds confirms Equation (66). Stars give 
the non-local contribution, appropriately normalized relative to the 
open diamonds, i.e., 5 NL = (H/4)((ri -r 2 )<5(«0 i i -«0,2)}c/(s g rid^) 2 
(Eqs. [49]— [50]). Bottom panel: post-collision relative velocities are 
indeed nearly isotropic, as we had surmised in Equation (63). For 
this figure, we take u = Vi — V2, the total relative velocity; it 
includes the difference in the circular Keplerian velocities at the 
locations of the two particles. 



these three terms can be measured from simulations. In 
making these and other measurements, we sought ways 
to minimize noise introduced by finite particle numbers 
(Poisson fluctuations). For example, when measuring vis- 
cous fluxes of angular momentum and energy, it proves 
useful to consider only those particles that actually collide 
during the measurement interval. 

The stage is now set for simulating more complicated 
systems — narrow eccentric rings (like the Maxwell and Ti- 
tan ringlets of Saturn, or the Epsilon ring of Uranus) , and 
circumstellar disks with embedded planets. Among the 
phenomena we are interested in exploring numerically are 
the formation of sharp edges by shepherd satellites, the 
evolution of narrow rings into states of rigid apsidal pre- 
cession, and the eccentricity evolution of planets as driven 
by disks. 

We thank Ruth Murray-Clay for helpful exploratory cal- 
culations, and Jack Wisdom for encouraging remarks. EC 
acknowledges support from the National Science Founda- 
tion, NASA, and the Alfred P. Sloan Foundation, and is 
grateful for the warm hospitality of the Canadian Insti- 
tute for Theoretical Astrophysics / University of Toronto, 
where a portion of this work was completed. 



7. SUMMARY AND OUTLOOK 

We have introduced an algorithm to simulate collisions 
between inelastic particles in an optically thin disk orbit- 
ing a central mass. The algorithm is simple to implement 
and adds negligible running time to existing collisionlcss 
N-body codes. A major feature of the algorithm is that the 
disk particles' motions need only be tracked in a plane. Yet 
the algorithm transcends its two-dimensional appearance 
to simulate a three-dimensional disk of particles whose 
random velocity distribution tends to be isotropized by 
collisions. 

We have performed a battery of tests of the algorithm 
for the case of an isolated, narrow, circular ring. Numer- 
ical simulations agree with analytic theory with regard to 
how the particles' velocity dispersion equilibrates, how the 
ring viscously spreads, how energy and angular momen- 
tum are transported, and how energy dissipation relates 
to the viscous angular momentum flux and to the back- 
ground shear. Angular momentum transport arises not 
only from particle advection (HF n ), but also from corre- 
lations in the random velocity field (F r( p) and from finite 
particle sizes (Fnl)- The relative magnitudes of each of 
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